Screening in 2D: GW calculations for surfaces and thin films using the repeated-slab 

approach 
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In the context of photoelectron spectroscopy, the GW approach has developed into the method 
of choice for computing excitation spectra of weakly correlated bulk systems and their surfaces. 
To employ the established computational schemes that have been developed for three-dimensional 
crystals, two-dimensional systems are typically treated in the repeated-slab approach. In this work 
we critically examine this approach and identify three important aspects for which the treatment 
of long-range screening in two dimensions differs from the bulk: (1) anisotropy of the macroscopic 
screening (2) k-point sampling parallel to the surface (3) periodic repetition and slab-slab interaction. 
For prototypical semiconductor (silicon) and ionic (NaCl) thin films we quantify the individual 
contributions of points (1) to (3) and develop robust and efficient correction schemes derived from 
the classic theory of dielectric screening. 
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I. INTRODUCTION 

Surface science has developed into a highly active and 
multidisciplinary area of research that has produced a 
rich variety of surface-sensitive experimental techniques. 
In this context theory and calculations have become an 
invaluable tool for interpreting the often complex and 
indirect data and for making predictions for structures 
or properties still inaccessible experimentally. 

For many experimentally and technologically relevant 
questions the electronic structure of a system is of cen- 
tral importance. Experimentally, it can be probed, for 
example, by photocmission spectroscopy. Theoretically, 
many-body perturbation theory in the GW approxima- 
tion, where G is the Green's function and W is the dy- 
namically screened Coulomb interaction, has developed 
into the method of choice for describing electronic exci- 
tations as measured in direct and inverse photocmission 
in weakly correlated solids and their surfaces.— i^>^>^>^>^>i 

Most GW implementations, notably such that employ 
pure or augmented plane- waves as basis functions f^i^iiSiiii 
but also others^ rely on three-dimensional periodic 
boundary conditions, which are appropriate for crys- 
talline bulk systems. To treat systems with a reduced 
periodicity like surfaces, thin films, nanowires, clusters, 
or molecules, three-dimensional periodicity is imposed ar- 
tificially. For this, the system of interest is placed in a 
three-dimensional unit cell, with an empty (vacuum) re- 
gion to separate the physical system from its periodic 
images in the broken-symmetry direction(s); see Fig. [TJ 
This procedure is frequently referred to as "supercell ap- 
proach," but for systems with a two-dimensional period- 
icity we prefer the more descriptive term "repeated-slab 
approach." 

In this paper we will focus on the two-dimensional 
case, i.e., the description of surfaces and thin films. Al- 
ready in the early days of modern GW calculations, Hy- 
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FIG. 1: Repeated-slab approach for surfaces and thin films 
schematically, with s being the slab thickness and v the slab 
separation. 



bertsen and Louie applied the then new methodology 
to simple semiconductor surfaces using the repeated-slab 
approachii^ In their paper they raised several fundamen- 
tal and technical questions. Of importance here is their 
speculation that "the crucial change in the self-energy 
operator at the surface may be largely contained in the 
Green's function," while "the screened interaction may 
more closely follow the variation in the local density." 
They further emphasize that their preliminary conclu- 
sions "will require further examples ... as well as a crit- 
ical evaluation of the slab approach for the surface self- 
energy operator" . Surprisingly, these issues have only 
rarely been addressed in later studies and the "critical 
evaluation of the slab approach" is still missing in the 
literature. This paper is a contribution to fill this gap, 
illustrating that the screened interaction does not simply 
follow the local density but unfortunately is substantially 
influenced by the repeated-slab geometry. 

In principle it is straightforward to investigate how 
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the artificial periodic repetition (in one or more dimen- 
sions) affects a given computational method by increas- 
ing the vacuum region until the properties of interest ex- 
hibit no dependence on this parameter anymore. This 
is easily achieved in density functional theory (DFT) 
with the commonly employed local or semilocal func- 
tionals, provided that appropriate correction schemes 
for the low-order electrostatic multipolc moments arc 
appliedfi^ if necessary. The decoupling becomes more 
difficult for GW due to the long-range nature of the 
(screened) Coulomb interaction, as has been explicitly 
shown for a Na4 cluster To solve the decoupling prob- 
lem, it has been proposed to cut off the Coulomb inter- 
action in the broken-symmetry directions . ^^i^^'^^i^^ How- 
ever, to ensure that no interactions within the physical 
system arc truncated, the vacuum region must be at least 
as large as the physical system in these schemes, which 
makes the approach unproportionally expensive (in terms 
of the computational effort) for thicker slabs. Moreover, 
the k-point sampling of the Brillouin zone, the compu- 
tational parameter naturally associated with these long- 
range effects, becomes a critical convergence parameter 
along the nontruncated directions F^^i^^ As we will show 
below, even this seemingly technical issue is related to 
the screening and hence the physical properties of the 
system. 

An obvious way to avoid these spurious interactions for 
systems with broken translational symmetry is to aban- 
don the concept of periodic boundary conditions alto- 
gether in the relevant directions and perform the calcu- 
lations entirely in real space. For semi-infinite jellium 
surfaces such a GW embedding scheme has been suc- 
cessfully implemented i^^'^° Its extension to realistic sur- 
faces, however, is computationally still too expensive. A 
real-space implementation for finite systems has also been 
report ed ,^^'^^ but its applicability to systems with peri- 
odicity in one or more directions remains to be shown. 

The purpose of this work is to reevaluate the perfor- 
mance of the repeated-slab approach when no Coulomb- 
truncation technique is employed. We address the nature 
and the magnitude of the effects introduced by the pe- 
riodic repetition and propose a robust correction scheme 
that allows to extract the isolated-slab limit already from 
very small vacuum separations. To achieve this, it proved 
necessary to explicitly take the anisotropy of the macro- 
scopic screening into account. In addition, wc have found 
that the k-point sampling of the Brillouin zone becomes 
a more critical parameter for GW repeated-slab calcu- 
lations than recognized so far. Appropriate samplings 
sufhcient for bulk GW or slab DFT calculations cannot 
be transferred to surfaces in general. 

For weakly correlated bulk systems GW calcula- 
tions are now routinely performed. Remaining open 
issues, such as computational efficiency, influence of 
pseudopotentials or self-consistency, are actively being 
addressed^iaa2i2i^i24^i26^ but do not affect the con- 
clusions drawn in this paper. For more strongly corre- 
lated systems, it becomes necessary to go beyond GW, 



but these schemes often include the GW self-energy dia- 
grams as lowest order<^2i2^ Similarly, the Bethe-Salpeter 
approach to electron-hole excitations, as probed in op- 
tical absorption or electron energy-loss spectroscopy, 
builds on the GW self-energy In these schemes, screen- 
ing plays a similar role as for GW calculations, but they 
go beyond the scope of the present work. 

We note also that the surface electronic structure of a 
slab may differ from that of a semi-infinite solid. Surface 
resonances and plasmons, for example, may be affected 
by confinement effects. The finite thickness of the slab 
becomes an additional potential source of error in sur- 
face calculations, but will not affect free-standing films 
or heterostructures. Converging the surface electronic 
structure with respect to slab thickness and/or develop- 
ing robust correction schemes is a separate problem that 
would require a detailed study in its own right. We will 
briefly address the relevant issues in Sec. IIIII 

The remainder of this paper is organized as follows: In 
Sec. In] we present the methodological and physical as- 
pects of the repeated-slab approach. After briefly sum- 
marizing the GW space-time method (Sec. Ill Ap . we 
focus flrst on the anisotropy of macroscopic screening 
(Sec. Ill Bp and the k-point convergence parallel to the 
surface (Sec. Ill Cp . Then we discuss the influence of the 
periodic repetition and the associated convergence of the 
band energies with respect to the vacuum separation be- 
tween the slabs (Sec.|TlD|). In Sec.lTlIl we summarize our 
main findings and put our conclusions into the context 
of other computational GW schemes. In the Appendix 
we present the computational scheme for the classical di- 
electric models employed. 



II. LONG-RANGE SCREENING 
A. The GW space-time method 

All calculations have been performed with the GW 
space-time method.— '^^i"^^ Due to its advantageous lin- 
ear scaling behavior with respect to the k-point sam- 
pling of the Brillouin zone it is ideally suited for per- 
forming the extensive convergence studies presented be- 
low. Wc have recently extended the code to include the 
anisotropy in the long-range screeningfS^, a crucial point 
for the repeated-slab approach, as we will demonstrate 
below. In the following, we will only briefiy sketch out 
the computational scheme and refer the interested reader 
to the previous papers for further details. 

In the space-time method, the Green's function G is 
constructed in real space (r and r') and imaginary time 
(ir) from the Kohn-Sham wave functions (/)„k and ener- 
gies e„k (the Fermi level is set as the energy zero). 

G(r,r';±zr) = ^/ d^fc ^ ^„k(r)<k(r')e-^-^ 
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where f2 denotes the unit-cell volume and the integral 
over k runs over the first Brillouin zone. Depending on 
the signs taken, the state summation over n runs over un- 
occupied (occupied) states for positive (negative) imagi- 
nary times. 

The polarizability is calculated in the random-phase 
approximation by 



P{r,r';iT) = -2iG(r, r'; iT)G'(r, r'; -ir) 



(2) 



and is then Fourier transformed to reciprocal space and 
imaginary frequency (iuj). The polarizability P is used to 
compute the screened interaction W via the symmetrized 
dielectric matrix 



ek(G,G';iw) = Sac 



An 



|k + G||k + G 



-Pk(G,G';*c.) 



(3) 

From its matrix inverse for each k and iuj, the screened 
interaction 



Wi,{G,G';iuj) 



Ik-hGllk + G 



-£-i(G,G';zc.) (4) 



is obtained and then Fourier transformed back to real 
space and imaginary time. The self-energy is then com- 
puted as 



E(r, r'; ir) = iG(r, r'; iT)W{v, r'-ir) 



(5) 



To obtain the self-energy corrections to the Kohn-Sham 
energies e„kj the matrix elements of the perturbation op- 
erator E — Vy^c (where T^c denotes the local exchange- 
correlation potential) are computed and Fourier trans- 
formed to the imaginary frequency axis. These matrix 
elements are then analytically continued to the real fre- 
quency axis by fitting a multi-pole function on the imag- 
inary frequency axis.— Finally, the quasiparticle energies 
^nk ^'^^ given by the solution of 



,qp 



Enk + ((/3„k|E(eJi^) - KclVnk) 



(6) 



where (iy9nk| ly'rik) denotes matrix elements with respect 
to the Kohn-Sham wave functions. 



B. Anisotropy of the macroscopic screening 

An important point for computing the screened inter- 
action in reciprocal space [Eq. ©-([I])] is the sin- 
gularity of the Coulomb potential as k ^ 0. For the 
G = G' = element of the symmetrized dielectric ma- 
trix [Eq. ^] it is cancelled by the behavior of the 
polarizability In practice, we enforce this cancellation 
analytically by performing a Taylor expansion for the 
polarizability around the F point k = 0. However, the 
expansion depends in general on the spatial direction in 
which the F point is approached. This directional depen- 
dence reflects the anisotropy in the macroscopic screening 
and introduces a nonanalytic, yet finite contribution to 
the inverse dielectric matrix 




FIG. 2: Dielectric tensor for repeated-slab systems according 
to effective-medium theory as a function of the s/c ratio (s 
slab thickness, c total size of the simulation box): component 
£|l parallel to the surface and Ez perpendicular to it. 



The screened interaction in the vicinity of the F point 
then takes the form^i^i^ 



Wk(0,0;^CJ) 



A-K 



kTL(icj)k ' 



(7) 



where L(icL') denotes the macroscopic dielectric tensor. 
The tensor expression in the denominator reduces to a 
scalar when the macroscopic screening is isotropic. To 
our knowledge, most GW implementations explicitly ex- 
ploit this simplification (exceptions being the exact treat- 
ment of Ref. |36| and the approximate but accurate^ ten- 
sor treatment of Ref. ISTt ) and therefore implicitly assume 
that the screening is isotropic. However, this is not the 
case for repeated-slab systems, as is easily demonstrated. 
Assuming that the slabs in Fig. [T] are homogeneous di- 
electrics of finite thickness s with an isotropic bulk di- 
electric constant eb and separated by a vacuum region 
with thickness w, the dielectric tensor components of the 
repeated-slab system are given by 



ehs + v s 
e\\ = ■ = 1 + (eb - 1)- , 

S + V c 



.-1 



•^ + ^'-l-(eb-l)^, 



S + V 



(8) 
(9) 



where c ^ s + v denotes the total height of the simulation 
cell. £|[ (parallel to the surface) and (perpendicular to 
it) agree for isotropic systems (s/c = or 1) but deviate 
considerably for s/c ratios between these limiting cases 
(cf. Fig. [21). The ratio e\\/ez becomes largest for s/c = 
1/2. In practice, the parameters of repeated-slab systems 
are often close to this ratio of maximum anisotropy. 

In the GW calculations the singular elements of W 
require a special treatment in the subsequent computa- 
tional steps, which formally involve an integration over 
the Brilloiun zone. Since the singularity is integrable, 
the result is formally well defined, but cannot be approx- 
imated by a finite summation. The solution is to split 
the screened interaction according to 



(10) 
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(a) isotropic 




FIG. 3: Convergence of the fundamental gap of a H-saturated 
four-layer Si(lOO) slab with respect to the number of k points 
A^z perpendicular to the surface for (a) the isotropic and (b) 
the anisotropic approach. The quantitative behavior depends 
on the sampling in the parallel direction (A'^y x Ai'n). Note the 
different scales of the two graphs. 



into a long-range part W^" and a short-range remainder 
^sr_ Pqj, ^ir^ a simple analytic form is chosen that ex- 
hibits the same singularity as W, but for which the inte- 
gral can be computed analyticallyiSii^i^i^i^^i^i^^ The re- 
mainder W^"^ then becomes nonsingular, and its integral 
can safely be replaced by the sum over a finite k-point 
grid. 

We demonstrate the importance of an anisotropic 
treatment for W^'^ for the case of a hydrogen-saturated 
four-layer Si(lOO) slab with a vacuum region equivalent 
to four layers of silicon. The calculated nonzero ele- 
ments of this repeated-slab system's dielectric tensor are 
£xx=5-^, Syy=5.5, and £2^=2.2 at the smallest imagi- 
nary frequency a;=0.036 hartree, highlighting the large 
anisotropy predicted by the electrostatic considerations 
above. In Fig. [3l we show the convergence of the fun- 
damental gap^ with respect to the k-point sampling 
perpendicular to the surface. It is obvious that isotropic 
averaging for the screened interaction, i.e., using a scalar 
rather than a tensorial expression for the singularity of 
W^'^ , leads to an unphysical linear increase in the band 
gap, which will not level off when is increased fur- 
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FIG. 4: Band gap for a two-layer NaCl(lOO) slab and a 
hydrogen-saturated four-layer Si(lOO) slab with a N\\ x A'^y x 1 
k-point sampling as a function of l/A'^y . The vacuum region is 
20 bohr wide for both systems. The fit function is explained 
in the text. 



ther. The reason for this linear divergence lies in the 
inadequate treatment of the singularity in Eq. ([7]) , which 
is not fully removed.^ In contrast, the anisotropic treat- 
ment converges rapidly. 

Only the proper anisotropic treatment allows us to in- 
vestigate the importance of the k-point sampling in the 
direction perpendicular to the surface. This convergence 
has not been discussed for repeated-slab systems before. 
Using the anisotropic treatment in the GW space-time 
method and a A^|| x A^|| x iV^ sampling, we find that the 
self-energy corrections exhibit a l/iV^ behavior, the mag- 
nitude of which rapidly decreases with increasing 
Such a behavior might result from the remaining approx- 
imations made for the F point. In practice, the conver- 
gence with respect to Nz plays a relevant role only for 
coarse ky samplings. We will address convergence with 
respect to this parallel sampling in Sec. Ill CI Test cal- 
culations indicate that the qualitative behavior is not af- 
fected by the choice of TV^ . For simplicity, we have there- 
fore used = 1 for the calculations reported in the 
following. 



C. k-point sampling parallel to the surface 

The k|| sampling is a further critical aspect in GW 
calculations for surfaces and thin films, which may not 
have attracted sufficient attention so far. We demon- 
strate the importance of this issue for two representative 
slab systems: a two-layer NaCl(lOO) slab and a hydrogen- 
saturated four-layer Si(lOO) slab. In Fig. [4]the gap for a 
A^ll X A^ii x 1 sampfing is plotted as a function of 1/^|| . In 
both cases we observe a similar shape of the curves that 
does not follow a simple power law as observed for iV^,. 
We will show in the following that this can be understood 
in terms of the particular decay behavior of the screened 
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FIG. 5: The nonlocality of two-point functions, such as E, is 
defined within the interaction cell. In the space-time method, 
its size is linked to the k-point discretization grid. 
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FIG. 6: a) Computation of the screened potential from the 
image charges, b) Periodic repetition of the image charges 
to simulate the effect of the ky discretization, a denotes the 
lateral extent of the interaction cell. 



interaction in a slab system and can be modeled using a 
simple analytic function. 

To better understand the approximations introduced 
by a finite k-point sampling, we will first derive the gen- 
eral connection between the k-point sampling and the 
range of nonlocality, and then discuss its relevance for 
the present case. As an introductory remark, we note 
that in periodic systems the two-point functions F = 
M^, or S reflect the lattice periodicity 

F(r,r') =F(r + R,r' + R) , (11) 

where R denotes a vector of the real-space lattice defined 
by the unit-cell vectors ai, a2, and aa. We can then 
introduce the representation 

FR(r,r') :=F(r + R,r') , (12) 

where r and r' are restricted to the unit cell. The cor- 
responding reciprocal-space representation Fk (also de- 
noted as mixed-space representation)^ is obtained from 
a Fourier transformation 

Fk(r, r') := e-^k.(r-r') ^ p^^^^ ^-.k-R (^3) 

R 

We denote the unit vectors of the reciprocal lattice by 
bi,? G {1,2,3}. They are defined by those of the real- 
space unit cell via a^bj = 27r5y. The connection to the 
R representation provides a real-space picture of the k- 
point discretization. For this, we note that a regular, 
F-centered discretization grid (iVi x N2 x N3) defines a 
lattice hi/Ni in reciprocal space that is associated with 
a real-space supercell NiRi, comprising Ni x N2 x N3 
unit cells (cf. Fig. [3]). This supercell coincides with the 
interaction cell of the space-time method;^ or in other 
words the range of nonlocality for the two-point function 
F = G, W, and S. Discretizing k for the function Fj. is 
therefore equivalent to imposing a translational symme- 
try 

FR(r,r') =FR+s(r,r') , (14) 

in real space, where S denotes a lattice vector of the 
interaction-cell lattice (iV^ai). 



From this connection we conclude that convergence in 
the k|| sampling can only be achieved when the inter- 
action cell associated with the sampling is large enough 
to encompass the dominant nonlocality in the relevant 
functions (G, P, W, and S). This also holds true for 
GW schemes that do not explicitly rely on the real-space 
interaction cell, or which employ a representation other 
than plane waves for the unit-cell coordinates (r, r'). Let 
us now consider the behavior of the different two-point 
functions. In nonmetallic bulk systems, the Green's func- 
tion decays exponentially with increasing distance of its 
arguments, typically over a few bond lengths4£ We see 
no reason to believe that repeated-slab systems show a 
qualitatively different behavior in this respect. However, 
the surface may exhibit a different band gap than the 
bulk, or it may even be metallic. The decay properties of 
the Green's function will then change accordingly .i^Sdi In 
particular, a reduction in the gap at the surface implies 
a slower decay of G at the surface compared to the bulk. 
The decay properties of G are directly transferred to P 
[Eq. (HD] and S [Eq. (O]. The screened interaction W, 
on the other hand, exhibits a very slow 1/r decay. How- 
ever, this long-range limit H^'"' is subtracted from the full 
W in the reciprocal-space singularity treatment discussed 
previously. For the numerical convergence only the re- 
mainder W^'^ is relevant, i.e., the difference between W 
and the model function W^'"' used in the singularity treat- 
ment. W^'^ is dominated by the variation of the electron 
density and thus by the atomic structure, and should 
become negligible when the typical length scale of these 
variations is exceeded. In bulk systems, this corresponds 
to a few bond lengths. However, the decisive structural 
variation in a repeated-slab system is the slab itself (at 
the slab boundary the dielectric constant drops from £b 
to 1). The thickness of the slab or the vacuum therefore 
defines the length scale for the interaction-cell conver- 
gence also parallel to the surface. 

To understand how the k|| discretization modifies the 
screened interaction, we employ the classical theory of di- 
electric screening. Similar to what was done above for the 
macroscopic dielectric tensor, we approximate the slab 
by a homogeneous dielectric with a sharp interface to 
the vacuum as sketched out in Fig. [HI We then apply 
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FIG. 7: (Color online) Additional potential AV™ introduced 
by the ky discretization of the screened interaction (see text). 
The dielectric model corresponds to a two-layer NaCl(lOO) 
slab in a c=30bohr unit cell with a 3x3 ky sampling. The 
blue-shaded areas indicate the position of the slabs; the hori- 
zontal dashed line marks the average potential. The isolated- 
slab case is shown for comparison. 

the method of image charges (cf. Appendix) to compute 
the screening function. When a charge is placed at point 
z' , the resulting image charges lie on a line through the 
original charge perpendicular to the surface (cf. Fig.lH]). 
The effect of the k|| discretization is simulated by impos- 
ing translational symmetry for the image charges in the 
direction parallel to the surface according to Eq. In 
other words, we explicitly consider a lattice of perpen- 
dicularly aligned image charges [cf. Fig. Eljb)] with the 
interaction cell's lattice constant of A^n times the original 
cell's lattice constant. To simplify the discussion we focus 
on one characteristic aspect of the screened interaction: 
the potential generated by the image charges at the 
position of the original charge, i.e., aX z = z' . We obtain 

where the sum over S runs over the lattice vectors of 
the interaction cell iVnan. d and a enumerate the image 
charges that arc computed as described in the Appendix. 
From this the long-range model function 1/ {e\\r) and the 
contribution of the original image charges at S = are 
subtracted. Exploiting the sum rule Eq. ([35]) . wc obtain 
for the error introduced by the parallel discretization 

In Fig. [7] we show AT^™(z) for the dielectric model 
corresponding to the two-layer NaCl(lOO) slab with a 3 x 
3 k|| sampling for both the isolated and the repeated- 
slab case. The absolute error in the potential is very 
large. We note that the average shift in the potential 
is implicitly corrected for in the full GW calculations by 
setting the VFo(0, 0) element to zero after the subtraction 



of the model function. Therefore, only the deviation from 
this average (dashed line in Fig. [7| contributes to the 
k-point convergence, which is considerably smaller but 
far from negligible. We conclude that the convergence 
of a GW calculation with respect to the k|| sampling 
is crucially influenced by the discretization error of the 
screened interaction. 

The variation of AX^™ from its average value may be 
taken as an indicator for the discretization error in the 
GW self-energy. We investigated how this changes as a 
function of slab and vacuum thickness within the dielec- 
tric model. When the vacuum region is increased, the 
variation increases because AF"" quickly decays in the 
vacuum region (cf. Fig. [71), thereby pulling down the 
average over the computational cell. The periodic rep- 
etition (small vacuum) thus introduces a highly advan- 
tageous compensation effect: While the absolute value 
of Ay™ is quite large in the slab, the average is domi- 
nated by the contribution of the slab, which then cancels 
out a large part of the total discretization error. The 
physical origin behind the slab and vacuum dependence 
of the k|| convergence lies in the scale dependence of the 
screening in such an inhomogeneous system. For dis- 
tances much smaller than the slab thickness, a bulk-like 
screening is expected inside the slab. However, in cells 
with large vacuum separations the components of the di- 
electric tensor may deviate considerably from this bulk 
limit; as illustrated in Sec. IIIBI Since the treatment for 
the long-ranged part of the screened interaction assumes 
this behavior at all length scales, a portion of this differ- 
ence in screening enters the numerical treatment of the 
short-ranged part W^^ and affects the k|| convergence. If, 
on the other hand, the vacuum separation is small, the 
dielectric screening of the repeated-slab system is close 
to bulklike and the magnitude of the discretization er- 
ror is reduced correspondingly. We find this anticipated 
behavior fully confirmed for the slab systems considered 
in this work (Figure [10] will show this explicitly for the 
two-layer NaCl(lOO) slab). 

To estimate the infiuencc of the W discretization on the 
quasiparticle energies computed in the GW space-time 
method more precisely, we have developed a fit function 
to describe the dependence of the band energies on the k- 
point sampling A^||. For this purpose, we assume that the 
Go Wo corrections show the same functional dependence 
on iV|| as AF™ in Eq. (fTB]) . Retaining a single term of 
the summation, we arrive at a three-parameter function: 

e„(Ar||) = .„(oo) + - • (17) 

The parameters Qn-, Dm and e„(oo) are determined for 
each state n by fitting to the GW data. We find that 
this relatively simple form accurately describes the con- 
vergence with respect to iVy for all slab systems studied. 
Typical examples of the quality of the fit are shown in 
Fig. [D We employ the fitting procedure to estimate the 
remaining error at finite k|| sampling, or to extrapolate 
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FIG. 8: Dependence of the band gap of H-saturated Si(lOO) 
slabs with four, six, and eight layers (two layers ~ 5 bohr) 
on the total cell height c. Solid lines: DFT+GW results. 
Dashed lines: DFT+GW results with finite-vacuum correc- 
tions derived from the dielectric model. The DFT gaps (not 
shown) are independent of the vacuum thickness. 



the converged value e„(oo) in cases where the error at 
finite samphng is unacceptably large. 



D. Periodic repetition: slab-slab interaction 

Until now, wc have focused the discussion on technical 
and numerical aspects of GW calculations for repeated- 
slab systems. We now turn to the physical effects asso- 
ciated with the interslab polarization. For this, we com- 
puted the band structures for hydrogen-saturated Si(lOO) 
slabs with four, six, and eight Si layers and varying 
amounts of vacuum. The k-point convergence was tested 
for each system, based on the fitting procedure described 
above. An 8x8x1 (7x7x1) sampling was found to 
be sufficient for the four-layer (six- and eight-layer) case. 
The numerical results for the fundamental gap presented 
in Fig. |8] show a slow convergence with respect to the 
vacuum size. The changes within the numerically acces- 
sible range are noticeable (ranging between 0.05 cV and 
0.2 eV per 10 bohr of vacuum). More importantly, wc es- 
timate that vacuum regions of 80-100 bohr would be re- 
quired to achieve absolute convergence within 0.05 eV of 
the isolated-slab limit (which can be determined with the 
correction scheme that will be described later). However, 
only the gap between occupied and unoccupied states is 
affected. The changes within the occupied or unoccupied 
part of the spectrum are minor. As we will show, this is a 
further effect of the long-range screening that affects the 
occupied and unoccupied bands as a whole: The pres- 
ence of the neighboring slabs gives rise to an additional 
contribution in the screened potential, which decreases 
when the distance between the slabs is increased. 

We investigated this within the dielectric model. The 
model parameters are obtained consistently with the full 
GW calculation from the dielectric tensor components e|| 
and Ez and the total cell height c. From Eqs. ([5]) and 
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(18) 
(19) 



To simulate the effect of the periodic repetition, we take 
into account a finite number of slabs («20) in our model 
calculation and compare this with the limiting case of an 
isolated slab, which corresponds to increasing the vacuum 
size to infinity. 

The change in the self-energy that results from the 
long-range screening effects can be estimated from 
the static Coulomb-hole screened-exchange (COHSEX) 
approximation^: 



S(r,r') = ScoH(r,r') + SsEx(r,r') , 



(20) 



EcoH(r,r') = -<5(r - r') ^r, r') - z;(r, r')] ,(21) 



(r,r') = -^^„(r)^;(r')M^(r,r') 



(22) 



Here W denotes the statically screened interaction and 
V the bare Coulomb interaction. To isolate the contribu- 
tion arising from the neighboring slabs we separate the 
screened interaction into a bulklike part T^^*^""^ and an 
additional contribution iyp°' from long-range polariza- 
tion effects. Since S is a direct product of W with the 
Green's frmction, the self-energy reflects this separation: 



occ 

I]P°'(r,r')= -<5(r-r')-^<^„(r)<(r') W^P°'(r,r'). 

n 

(23) 

To simplify this expression further, we split off purely 
local contributions from W^°^ according to 

WP°\r,r') = ^[WP°\r,r) + WP°\r',r')] (24) 



W^P°''"'(r,r') 



(25) 



implicitly defining the purely nonlocal remainder VI^p°i^"'. 
We now argue that the self-energy 5]p°''"' arising from 
this separation is negligible for long-range polarization 
effects. In particular, we have by construction. 



iyP°''"'(r,r) = 



(26) 



The corresponding COH contribution therefore vanishes. 

Assuming that W^°^ is generated 



We now turn to S 



pol.nl 
SEX 



by a set of image charges at distance d, WP°^''^^{r, r') w 
for |r — r'l ^ c?. The magnitude of d is determined by the 
distance to the relevant dielectric interface. On the other 
hand, the second factor in the SEX term is the Green's 
function at t — > 0, t > 0. Its exponential decay limits the 
magnitude of |r — r'| to a few interatomic distances4S 



Therefore. S 



poLnl 
SEX 



is expected to be small for polariza- 



tion effects taking place at a length scale of more than 
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FIG. 9: Image potential for an isolated slab and a repeated- 
slab system in comparison. The blue-shaded regions indicate 
the position of the slabs. The model parameters correspond 
to a two-layer NaCl slab. 
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FIG. 10: Left: ky convergence for the gap of a two-layer 
NaCl slab for different total ceU heights c. Right: DFT-LDA 
gap (bottom) and extrapolated quasiparticle correction Aqp 
as well as its value corrected for the finite vacuum (top) as a 
function of c. 



a few bohr, and we can safely neglect its contribution. 
For SPoi^iot: arising from the local part [Eq. ([HI)], the 
state summation reduces to a projection onto the occu- 
pied states, and the expectation values of spoi^iot^ become 

(^„|SP°''i-|^„) ± i(^„|M^P°i(r,r)|^„) , (27) 

where the plus (minus) sign applies to unoccupied (oc- 
cupied) states. For the dielectric model, we now identify 
M^^""^withl/(£(z)|r-r'|). W^P°i(r,r) then reduces to the 
potential y™(z) induced by a charge at its own position, 
the S = term in Eq. (fT5|) . Equation (P7)) corresponds 
to an adiabatic switching on of "image charge" effects for 
the charged final state, taking into account the sign con- 
vention of single-particle energies. We note that Dclcruc 
et al^ have derived the same final formula from an elec- 
trostatic model and used it to estimate self-energy shifts 
in isolated nanoparticles. 

In Fig.m we compare T^™(z) for an isolated slab with 
that of the repeated slabs. We find for both cases a posi- 
tive contribution inside the slab, and a negative one out- 
side, in agreement with Delerue4^ The divergence at the 
interface is an artifact of the steplike dielectric profile em- 
ployed. In the repeated-slab case, however, the potential 
is shifted downward because of the additional polariza- 
tion in the neighboring slabs. It is precisely this polariza- 
tion of the neighboring slabs that creates the undesired 
perturbation in the central slab and is responsible for 
the observed dependence of the band gap on the vacuum 
size. The shift in the image potential (indicated by the 
dashed line in Fig. [5]) is essentially constant over the slab 
and continuous across the slab-vacuum interface. The 
EP°' self-energy corrections according to Eq. p7|) there- 
fore do not vary significantly between different states, 
giving rise to a scissor-like change in the gap in agree- 
ment with the observations made for the full GqWq cal- 
culation. Moreover, this fortuitously simple situation fa- 
cilitates a quantitative comparison between the dielec- 
tric model and the dependence of the band gap found 



in the GqWq calculation. To demonstrate this, we have 
extracted the magnitude of the finite- vacuum effect from 
the dielectric model and corrected the numerical GqWq 
data by these values for every vacuum thickness. The 
result has been included in Fig. [5] as the dashed curves. 
We find that this corrected data no longer depends signif- 
icantly on the vacuum size. Using this correction scheme, 
it becomes possible to determine the isolated-slab values 
for arbitrary vacuum sizes. 

These a posteriori corrections do not depend sensi- 
tively on the microscopic details of the material. We 
demonstrate this for an ionic material, a two-layer 
NaCl(lOO) slab. The dependence of the k|| convergence 
behavior on the slab/ vacuum ratio (left side of Fig. [T0|) 
agrees very well with the predictions from the dielectric 
model: With increasing vacuum size, the convergence 
becomes more difficult. For this reason, we found it nec- 
essary to extrapolate the ky convergence for each vacuum 
thickness using Eq. pT|) , because its magnitude becomes 
comparable to that of the finite- vacuum effect. Using 
the extrapolated values, we find a very good agreement 
for the finite-vacuum effect compared to the dielectric 
model. This is demonstrated by the quasiparticle correc- 
tion Aqp to the band gap shown at the right-hand side 
of Fig. 1101 Including the finite-vacuum correction, Aqp 
coincidcntally converges even faster with respect to the 
vacuum size than the band gap of the underlying DFT- 
LDA calculation in this case. 



III. DISCUSSION AND CONCLUSIONS 

In this work we have investigated the performance of 
the repeated-slab approach in GqWq calculations. Using 
the classic theory of dielectric screening we have found 
that the relevant effects can be understood from sim- 
ple dielectric slab models. We have demonstrated that 
long-range polarization effects introduce several impor- 
tant differences in slab systems compared to the cor- 
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responding bulk systems for prototypical semiconductor 
(Si) and ionic (NaCl) materials. The periodic repetition 
of the slabs introduces an additional scissorlike change 
of the one-particle spectrum that can be quantitatively 
corrected for a posteriori. These corrections are derived 
from dielectric slab models. 

Since long-range effects are described by the small-k 
behavior in reciprocal space, this limit requires special 
attention for slab calculations. The anisotropy of the 
macroscopic screening at k ^ and the slow conver- 
gence with respect to the k|| grid must be taken into 
account. If this is not done correctly, it may not be pos- 
sible to assess the importance of the vacuum separation. 
However, a careful treatment of these aspects does not 
pose principal problems and can be achieved with moder- 
ate computational effort, in particular for small vacuum 
sizes. Using the finite- vacuum correction scheme provides 
a highly efficient technique to extract accurate values for 
the isolated-slab limit from only a few calculations for 
small vacuum separations. 

We will briefly discuss how far the issues raised in this 
work for the GW space-time method are relevant also 
for other computational schemes. The anisotropy of the 
dielectric tensor is specific to the repeated-slab approach 
when no Coulomb-truncation techniques are employed. 
We have demonstrated recently^ how the anisotropy 
can be treated to arbitrary precision in plane-wave ap- 
proaches by expanding the angular dependence in spher- 
ical harmonics. In practice, already a maximum angu- 
lar momentum of 2 is sufficient to capture most of the 
effect 1^ We expect that a similar strategy is applicable 
for other basis sets, too. 

The slow convergence with respect to the k|| sam- 
pling is rooted in the reciprocal-space computation of 
the screened interaction. To our knowledge, all GW im- 
plementations with periodic boundary conditions rely on 
this strategy. However, studying the convergence with re- 
spect to this parameter is computationally very efficient 
in the space-time approach since it scales linearly in the 
number of k points, compared to a typically quadratic 
scaling in convolution approaches. The discretization er- 
ror arises from the part that is treated numerically, i.e., 
the deviation of the screened interaction from the model 
function used for the 1/fc^ singularity treatment. Al- 
ternative integration schemes, such as the improved in- 
tegration scheme by Pulci et alM- or the offset F-point 
methodfi^ will probably not improve this significantly, 
because they aim at better integrating the l/(k'^Lk) 
model function, not the deviations from it. The k|| con- 
vergence is a critical issue also for isolated slabs when 
treated by Coulomb-truncation techniques. In fact, it 
becomes even more dramatic in these cases because the 
long-range tail then equals the unscreened 1/r Coulomb 
interaction and the small-vacuum compensation effect is 
excluded. Indeed, Coulomb-truncation techniques have 
been found to require drastically enlarged k-point grids 
in the nontruncated directions J^J^ii^ We note that the 
anisotropy and the k-point convergence are relevant also 



for condensed slablike systems, such as multilayers or 
quantum- well superlattices, but their magnitude depends 
critically on the dielectric constants of the materials in- 
volved. 

The dependence of the gap on the vacuum thickness, 
on the other hand, is a physical property of repeated- 
slab systems. Since not only the closest neighboring 
slabs contribute to this effect, its magnitude reduces with 
increasing slab thickness4^ It is reproduced with the 
model calculation and is thus independent of the com- 
putational scheme used for the GqWq calculations. Any 
implementation that does not find this trend must em- 
ploy additional approximations that suppress this behav- 
ior. Whether such an implementation then provides the 
isolated-slab limit depends on the actual approximations 
and cannot be foreseen in general. 

For future GW calculations in repeated-slab systems, 
we recommend the following procedure: 

1. Use an anisotropic treatment of the screened 
Coulomb singularity for systematic k-point conver- 
gence studies. 

2. Check the k-point sampling parallel to the surface 
separately for each slab and vacuum thickness. If 
necessary, extrapolate using a physically justified 
function. 

3. Correct band energies/gaps for the finite- vacuum 
thickness according to Eq. ((27|) and the Ap- 
pendix. Equations (|18p and ([T9|l provide model pa- 
rameters consistent with the actual dielectric ten- 
sor. 

For surface calculations, it should be kept in mind that 
the surface electronic structure of a slab differs from that 
of a semi-infinite solid due to confinement of the elec- 
tronic states. This aspect is already present at the level 
of DFT and is directly transferred to the GW band struc- 
ture. Quantum confinement may for instance split a sur- 
face resonance into a series of quantized states. In addi- 
tion states localized at the two slab surfaces may couple 
through the slab at insufficient slab thicknesses. GW 
corrections might become important if they shift the en- 
ergies of surface states or resonances relative to the bulk, 
as this alters the decay behavior. 

For the surface GW self-energy itself, the role of the 
slab approximation has, to our knowledge, not been stud- 
ied. An in-depth analysis of this question goes beyond 
the scope of the present paper, but we will briefly dis- 
cuss the relevant aspects for G and W. For semiconduc- 
tors the exponential spatial decay of the Green's func- 
tion leads to a rapidly decreasing influence of the second 
surface with increasing slab thickness. However, if this 
decrease is fast enough to be unimportant at the com- 
monly employed slab thicknesses of only a few atomic 
layers cannot be guaranteed a priori. For the screened 
interaction, on the other hand, the influence of the di- 
electric discontinuity decays only slowly. Using simple 
electrostatic models similar to those detailed in Sec. Ill CI 
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(see also Delerue et alM-), we find that the presence of 
the additional surface at a distance s gives rise to an 
image-charge- like contribution w 2sh{eb+i)s ' Unlike for 
the finite- vacuum effect (Sec. IIID|) . however, quantita- 
tive corrections for this finite-slab screening effect depend 
strongly on the spatial extent of the electronic states. 

In summary, we find that the rcpcatcd-slab approach 
provides a computationally efficient way to calculate elec- 
tron and hole energies in the GqWo approximation for 
two-dimensional systems like surfaces and thin films. 
Three manifestations of long-range screening effects in 
two dimensions — the anisotropy of the macroscopic 
screening, the k|| sampling and the periodic repetition of 
the slabs — have been identified and analyzed in terms 
of the classic theory of dielectric screening. Their ef- 
fect on Go Wo calculations in the repeated-slab approach 
has been demonstrated for thin Si and NaCl films. Ro- 
bust and efficient correction schemes have been devel- 
oped. With these, the isolated-slab limit can be easily 
obtained from repeated-slab calculations. 
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COMPUTATION OF IMAGE CHARGES IN 
MULTILAYER SYSTEMS 

In this section we derive a practical scheme to compute 
the screened interaction in a dielectric multilayer system 
using image charges (cf. Fig. [6]). The screened interac- 
tion W{r, r') is obtained as the additional potential V{r) 
induced when a unit charge is placed at r'. V{r) is con- 
structed with a series of image charges vertically aligned 
above and below the original position. 

Let us first consider the textbook situation of two semi- 
infinite dielectric media Qi and Q2 with dielectric con- 
stants £1 and £2'^ For a charge q at r' in ili , the potential 
V{r) is given by 



1 



rel^i: y(r) = — 
r G : V"(r) 



£1 \ |r — r' I |r — r" 

1 q' 
£2 |r - r'l 



(28) 
(29) 



where q' and q" arc image charges, r" is obtained by 
refiecting r' at the interface, and we will therefore de- 
note q" as "refiected charge," whereas the effect of q is 



propagated into by the "propagated charge" q' . They 
are determined from the continuity equations of the elec- 
tric field and the electric displacement at the interface, 
yielding 



2£9 



q = 



£1 +£2 



and 



„ £1 — £2 

q = — : — q 

£1 +£2 



(30) 



In order to develop a computational scheme for a multi- 
layer system, a proper book-keeping is crucial to track the 
refiection and propagation of the various image charges. 
For this purpose, we divide the system into individual 
layers. Each layer has the same thickness L and a layer- 
specific dielectric constant £(z).— For the model calcu- 
lations in this work, we use L=l bohr. We will restrict 
the notation to the z coordinate for the position of the 
charges and layers. The parallel coordinate p becomes 
only relevant for the computation of the potential; see 
Fig. [nija). The origin of the coordinate system is chosen 
such that the layers are centered around integers L, and 
the interfaces arc at half integers. 

For reasons that will become clear below, we denote an 
image charge that contributes to the potential in layer 
z and is located at z + ad by q{z,d,a), where d > 
is the distance from the layer and cr = ±1. To show 
that the image charges can be determined iteratively we 
will now derive the iteration for d —>■ d + L. Consider 
a charge q{z,d,<7) relevant for the potential in layer z. 
Due to the interface at z — ^'^L two additional image 
charges appear. The reflected image charge, located at 
z — a(d + L), describes the potential in layer z and is 
given by [cf. (Eq. ^] 

q^'{z, d + £, -a) = ^i^j ~ ~ q{z, d, a) . (31) 
e{z) + e[z — a) 

The propagated image charge remains at the position 
z + ad and describes the potential in layer z — crL. Using 
our book-keeping notation and Eq. ([30|) , it can be written 

as 

q^^{z-a,d + L,a)^—^^^^^^q{z,d,a). (32) 

Obviously, the distance parameter is increased by L for 
each interface taken into account. The image charges 
for the distance d + L can thus be computed iteratively 
from those at distance d and will in general combine a 
reflected and a propagated contribution q ^ q"^ + qP''. 
The iterations are started by setting 



(Z(z',0,±l) = 1 



(33) 



We restrict z' to integer L, i.e., the original charge is 
placed at the center of a layer, thereby avoiding the di- 
vergence of the image potential at the dielectric discon- 
tinuities between adjacent layers. In practice, the itera- 
tions are stopped at some dmaxj which thereby becomes 
a convergence parameter. In addition, we truncate the 
number of layers considered explicitly and neglect image 
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charges that fall outside. This truncation becomes a sec- 
ond convergence parameter. The convergence for both 
parameters was tested by doubling the parameter until 
the changes became negligible. 

Summing the Coulomb potential of all the image 
charges relevant for this layer, one obtains 

The first term in Eq. p4|) describes the potential of the 
original unit charge for z = z' and corresponds to that 
of a bulk material with the local dielectric constant £{z). 
The sum over image charges q{z,d,<7) then introduces 



the long-range polarization effects due to the variations 
in the dielectric constant. 

As a final remark, we mention a useful connection to 
the macroscopic anisotropy. The long-range limit in the 
direction parallel to the surface l/(e|| p) is obtained for 
p ^ d by neglecting the perpendicular distance d in the 
denominator of Eq. p4p . From this, we obtain the sum 
rule 

-i^J*„,+|:,(..i.,j=i. (35) 
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